Read in data

Subsets only the NK and T cells.

library(Seurat)
## Attaching SeuratObject
## Attaching sp
library(dittoSeq)
## Loading required package: ggplot2
library(readxl)
sc <- readRDS("analysis/sc_integrated_annotated.rds")
tc <- subset(sc, subset = celltype_fine %in% c("NK cells", "T cells"))

rm(sc)

# load in file that can be used to re-orde samples in e.g. barplot
sample_order <- read_excel("metadata/sample_order.xlsx")
sample_order$sample <- sapply(sample_order$`Sample ID`, function(x){
  strsplit(x, " ")[[1]][1]
})

Re-do dimensiontality reduction on cycling cells only

DefaultAssay(tc) <- 'integrated'
tc <- ScaleData(tc, verbose = FALSE)
tc <- RunPCA(tc, verbose = FALSE)
tc <- RunUMAP(tc, dims = 1:50)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session

UMAPS

library(dittoSeq)
dittoDimPlot(tc, var = 'orig.ident', size = 0.25)

dittoDimPlot(tc, var = "Phase", size = 0.25)

dittoDimPlot(tc, var = "celltype_fine", size = 0.25)

Clustering

tc <- Seurat::FindNeighbors(tc, dims = 1:50, reduction = "pca")
tc <- Seurat::FindClusters(tc, resolution = seq(0.5, 1, 0.1))
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 35543
## Number of edges: 1811733
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8738
## Number of communities: 23
## Elapsed time: 12 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 35543
## Number of edges: 1811733
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8641
## Number of communities: 23
## Elapsed time: 12 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 35543
## Number of edges: 1811733
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8552
## Number of communities: 25
## Elapsed time: 12 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 35543
## Number of edges: 1811733
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8473
## Number of communities: 26
## Elapsed time: 12 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 35543
## Number of edges: 1811733
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8401
## Number of communities: 26
## Elapsed time: 12 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 35543
## Number of edges: 1811733
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8324
## Number of communities: 26
## Elapsed time: 13 seconds
pdf("plots/umap_t_cells_cl0.8.pdf", width = 9)
p <- dittoDimPlot(tc, var = 'integrated_snn_res.0.8', size = 0.5, do.label = TRUE)
print(p)
dev.off()
## quartz_off_screen 
##                 2
print(p)

tc$tcell_annotation <- "other"
DefaultAssay(tc) <- "RNA"

Plot marker genes

markers <- read.csv("metadata/marker_genes_man_tcell.csv")
mli <- split(markers, markers$tcell_annotation)

for(cellt in names(mli)) {
  marker_genes <- mli[[cellt]]$gene
  cellt_name <- gsub("[ \\+]", "", cellt)
  
  pdf(paste0("plots/UMAP_tcells_markers_", cellt_name, ".pdf"))
  multi_dittoDimPlot(tc, var = marker_genes,
                   min.color = "grey",
                   max.color = "purple",
                   size = 0.25,
                   order = "increasing")
  dev.off()
  
}
dittoDotPlot(tc, vars = unique(markers$gene), group.by = "integrated_snn_res.0.8")

Assign clusters to T cell type

tc$tcell_annotation[tc$integrated_snn_res.0.8 %in% c(3, 12)] <- "NK cells"
tc$tcell_annotation[tc$integrated_snn_res.0.8 %in% 10] <- "T regulatory cells"
tc$tcell_annotation[tc$integrated_snn_res.0.8 %in% 9] <- "T effector cells"
tc$tcell_annotation[tc$integrated_snn_res.0.8 %in% c(5,8)] <- "Naive T cells"
tc$tcell_annotation[tc$integrated_snn_res.0.8 %in% c(1,7,15,4)] <- "Cytotoxic T cells CD8+"
tc$tcell_annotation[tc$integrated_snn_res.0.8 %in% c(2,13)] <- "Exhausted T cells CD8+"

Visualization annotation

p <- dittoDimPlot(tc, var = "tcell_annotation", do.label = TRUE, size = 0.25,
                  labels.size = 2)
print(p)

pdf("plots/umap_manual_annotation_tcells_labeled.pdf", width = 10)
print(p)
dev.off()
## quartz_off_screen 
##                 2
p <- dittoDimPlot(tc, var = "tcell_annotation", size = 0.25)
print(p)

pdf("plots/umap_manual_annotation_tcells.pdf", width = 10)
print(p)
dev.off()
## quartz_off_screen 
##                 2
p <- dittoBarPlot(tc, var = "tcell_annotation", group.by = "orig.ident",
                  x.reorder =   match(sample_order$sample,
                                      unique(tc$orig.ident)))
print(p)

pdf("plots/barplot_manual_annotation_tcells.pdf", width = 10)
print(p)
dev.off()
## quartz_off_screen 
##                 2
library(dittoSeq)

p <- dittoDotPlot(tc, vars = unique(markers$gene),
                  group.by = "tcell_annotation",
                  min.color = "lightgrey", max.color = "darkred")
print(p)

pdf("plots/dotplot_tcell_markers.pdf")
print(p)
dev.off()
## quartz_off_screen 
##                 2
tmark <- read.delim("metadata/tcell_markers.txt", header = F)

cr <- colorRampPalette(c("white", "black", "red"))

p <- dittoHeatmap(tc, unique(c(markers$gene, tmark$V1)), 
                  scaled.to.max = TRUE,
                  heatmap.colors.max.scaled = cr(25),
                  annot.by = c("tcell_annotation", "integrated_snn_res.0.8"))

pdf("plots/heatmap_tcell_markers.pdf", height = 10, width = 15)
print(p)
dev.off()
## quartz_off_screen 
##                 2